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Abstract 

We compute isotope independent first and second order corrections to the Born-Oppcnheimer 
approximation for water and use them to predict isotopic shifts. For the diagonal collection, we 
use icMRCI wavefunctions and derivatives with respect to mass dependent internal coordinates 
to generate the mass independent correction functions. For the non-adiabatic correction, we use 
scaled SCF/CIS wave functions and a generalization of the Handy method to obtain mass indepen- 
dent correction functions. We find that including the non-adiabatic conection gives significantly 
improved results compared to just including the diagonal correction when the Born-Oppenheiiner 
potential energy surface is optimized for H 2 16 0. The agreement with expei imental results for deu- 
terium and tritium containing isotopes is nearly as good as our best empirical correction, however 
the present correction is expected to be more reliable for higher, unchaiacteiized, levels. 
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I. INTRODUCTION 


Starting with the work of Polyansky et al.,[ 1] empirical potential energy surfaces of water 
had reached a level of accuracy which enabled one to detect the need for mass dependent 


corrections arising from the breakdown of the Born-Oppenheimcr approximation. The first 
attempt to predict the mass corrections was carried out by Zobov et al.[ 2], who calculated 
the diagonal (first order) correction at the SCF level. They used this with an ab initio PES 
and found some improvement, although the errors were still too large to be unambiguous. 
Subsequently we[3] determined an empirical PES which reproduced H - 2 O energv levels 
very well, but when we applied the Zobov et al. correction, wn found that it did not predict 
the necessary isotopic shifts. Thus we determined an empirical mass dependent correction 


to our PES using data from HDO and D 2 0, and this greatly improved the results for all 
the other isotopoiners, however there remained discrepancies for the tritium isotopes, and 
because of the paucity of data for D 2 0, this correction is only reliable for low lying levels. [3] 


It would be most advantageous to the analysis and assignment of experimental spectra 


if reliable ab initio predictions of the mass dependent corrections is made available. This 
is because to a large extent, all of the information about the Born-Oppenheimcr PES is 


contained in the results of a single isotopomer. Experimentally, it is much simplei to deal 


with a single isotopomer. For example, one can obtain essentially pure samples of H 2 0 
and just obtain the spectrum of that species. However, it is not possible to obtain a puie 
sample of HDO, for collisions generate H 2 0 and D 2 0 as well, and so one has to deal with 
interfering H 2 0 and D 2 0 lines on top of the more complex spectrum that HDO has due 
to its lower symmetry. Molecules such as O 3 give rise to even more complicated mixtures 
of isotopoiners. Analyzing and reliably assigning experimental spectra is a great deal of 
work. If theoretitions could reliably predict isotopic shifts, then that would be a muc h more 


efficient means of generating spectral data. 

We compute the mass dependent corrections to the Born-Oppenheimer approximation 
using the second order perturbation theory result of Bunker and Moss [4] with ab initio wave- 
functions to compute the required matrix elements. [5] The first order correction is a mass 
dependent correction to the PES, sometimes called the Born-Oppenheimer diagonal correc- 
tion (BODC). The second order correction leads primarily to a mass dependent correction 
to the nuclear kinetic energy operator. This is usually called the noil-adiabatic correction. 
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We have recently developed tools to compute these corrections and applied them to H 2 0.[5] 
There we found that the SCF approximation used by Zobov to compute the first order 
correction is not reliable. In the present work we report predictions of isotopic shifts for 
deuterium and tritium containing isotopes using a newly optimized Born-Oppcnheimer FES 
including non-adiabatic effects and a more accurate first order correction. We also deiive a 
generalization of the Handy method for the diagonal correction to compute non-adiabatic 
correction functions from mass independent cartesian derivati ve matrix elements. 

II. OPTIMIZING THE PES 

In our previous work, [3] we obtained the empirical PES V' mp as 

yemp yZ y5Z _|_ ^core ^ycore _|_ ( .basis ^ybasis ^ ^\ rrvsl , ( 1 ) 

with c 5Z , c C0re , c basts , and the 18 parameters in AV /rest determined by a. least squares fit 
to 5493 ro-vibrational transitions with J < 5. This calculation did not include any non- 
adiabatic effects. In the kinetic energy operator, we used the nuclear masses. For the first, 
step in the present work, we re-optimize C' z , c rore , c basts , and the 18 parameters in AV 
including non-adiabatic effects. This produces the PES we call V nemp . 

The non-adiabatic corrections were obtained as an extension of our previous work. [5] 
Since the functional form used for the corrections does not extrapolate properly and in the 
present work we are considering higher energies, we extended cur grid in R\ and Ri to 0.775- 
1.425 A in steps of 0.05 A, and 0 to 50 - 170 degrees in steps of 10 degrees, and the energy 
cut off was 20,000 cm -1 . We fit these (693 unique) 1292 points to the same polynomial as 
in our previous work. In these calculations, the ground electronic state is treated using the 
SCF method, the excited states by the singles Cl method, and the results were scaled by 
1.1. This scale factor was obtained by matching the computed and experimental rotational 
g factors. [5] We used the aug-cc-pVTZ basis. [6] 

In our optimizations of the parameters in our PES, we use different experimental data 
than we used previously. In previous work we fit ro-vibrational transition eneigies. One 
problem with this is since ro-vibrational transition energies are differences in energy levels, it 
is hard to judge the amount of unique data one is fitting. Furthermore, there is the possibility 
of inconsistencies in ro-vibrational energies that make up the transition energies. So in the 
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present work, we use a different strategy. We start by compiling a list of ro- vibrational 
energies. These were determined as follows. We took the transition frequencies given in 
HITRAN 96, [7] and augmented this with the more recent data from the literature. [8 -10] 
We then iteratively solved for ro-vibrational energies with respect to the zero point energy. 
This gave many estimates of the ro-vibrational energies, and we only retained consistent, 
energies, which were defined as energies with at least 2 determinations, and a root-mean- 
square deviation from the average of the determinations less than 0.001 cm *. 1 his resulted 
in a total of 1688 energy levels having J in the range 0-8. In future work, we would use 
the compilation of Tennyson et al. [11] Now one desirable feature of our previous fit was that 
the errors had very weak dependence on J. To retain this feature, we weight rotational 
energies more than vibrational energies. We achieve this by organizing the ro-vibrational 
levels into groups which all have the same vibrational quantum number assignment. Then 
in each group, we subtracted the J = 0 energy from the higher J energies. If the .7 = 0 
energy has not been determined, then the lowest J level that has been determined is used 
instead. In the fit then, there is a relatively small number of “purely" vibrational energies 
and a much large number of rotational energies. We fit these with equal weights. 

For the present work we use a subset of the above data. This is to reduce the cost 
of the calculations and because we are not using hyperspherical coordinates. [3] We find 
our algorithm has difficulty converging highly excited states when we use r [ and 7 2 as 
coordinates. [5] Thus we consider only levels which correspond to roots 1-30 in a given JPS 
block. This results in a total of 1060 levels, of which 19 are “purely” vibrational levels and 
1041 are rotational levels. Wc made our best attempt to ensure the proper matching of 
experimental and theoretical levels, but inadvertently some miss-assignments were made. 
We attempted to automatically detect these during the least squares process by computing 
the average absolute deviation from the average of the difference between the observed 
and calculated levels. We multiplied this by ten, and used this as a cut-off criterian to 
greatly diminish the weights for levels which had greater difference between the observed 
and calculated levels. In our final fit, 33 points were deleted by this rule, and the weighted 
root-mean-square (rms) error of the fit was 0.054 cnr 1 . In this fit, we used the PES V emp 
and the kinetic energy operator including noil-adiabatic effects to generate the wavefunctions 
used as basis functions in the least squares fit. 

In the course of this work, we also carried out optimizations not including non-adiabatic 
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effects for comparison. We found only slight changes to our previous PES,[3] and that the 
errors produced with and without non-adiabatic corrections were essentially the same. The 
only differences were in the values of the parameters. Of great interest was that when we 
included non-adiabatic effects, the optimum value of the parameters c bai,s and c C0re came out 
to be essentially unity, and so in the final fit, we constrained them to be unity. It should be 
noted that we did not include the BODC or relativistic corrections! 12] to the PES in these 
optimizations. Initial trials including these corrections yielded essentially zero coefficients 
for the corrections. However, those calculations were carried out using the Zobov et al.[ 2] 
BODC, so if the calculations were repeated with the accurate BODC, this might not happen. 

III. BODC 

In Ref. 5, we reported calculations of the BODC of water using both CASSCF and SCF 
wavefunctions. For bending levels, the CASSCF and SCF results were very similar, but 
for the stretching levels, there were big differences. In the present work we have extended 
our capabilities and report the BODC computed at the icMllCI level, which should provide 
very accurate results. See Ref. 5 for the equations involved in the calculations. These 
calculations are carried out with a modified version of MOLPR02000.1.[13] Two additional 
modifications were required to carry out these calculations. First of all. we had to extend the 
diabatie orbital scheme to perform rotations among the virtual orbitals. Secondly, we had to 
evaluate two-particle expectation values using the icMRCI wavefunction. This was achieved 
by modifying the Cl code to perform a Cl iteration using the converged wavefunction, but 
using different one and two-electron integrals. This is complicated by the fact that the 
two-electron repulsion integrals are symmetric with respect to interchanging the indices 
for the basis functions of one electron, while for the two-particle operators we require, the 
integrals are skew symmetric with respect to interchanging the indices. Nonetheless we have 
successfully made these changes. 

When computing the BODC, one has three choices. Oik: can use the Handy method, 
[14] which for a triatomic requires nine derivatives, six of which can be evaluated using 
C s symmetry, and three which must be run in C\ symmetry. Alternatively, one can carry 
out calculations using internal coordinates that depend on the isotopomer. For example, 
for Radau coordinates, this requires three derivatives and all the calculations can be run 
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using C s symmetry. [5] Then the cost per point is less, but one has to perform additional 
calculations for each isotopomer of interest. Since in the present work we are interested in a 
large number of isotopomers, this would seem to favor the Handy method. However, avoiding 
the Ci calculations can be a big advantage. We therefore seek economies by combining the 
methods. This is done as follows. We know that for isotopomer ABC. the Handy method 
gives the BODC as 

AV = — F A H Fr H Fc, (2) 

m A m B rn c 

where m A is the mass of nuclei A, and F A is the correction function of atom A, etc. If 
we knew the BODC for three linearly independent isotopomers, we can solve for the F a , 
where a is A, D, or C. Specifically, we will carry out calcuations of the BODC using Radau 
coordinates for H 2 16 0, H 2 18 0, and HD 16 0. For the first two isotopomers, we only carry 
out calculations with the first HO bond length greater than the second, while for the last 
isotopomer, we do not make this restriction. We then solve for the F a at each geometry, 
and fit them to the functional form 

F a = exp{-/?[(r 1 -r e ) 2 + (r 2 -re) 2 ]}^c^[(r 1 -r e )/r e ] i [(r 2 -r e )/r c ] j [cos(0)-cos(0 e )] /c (3) 

ijk 

and then use Eq.2 to evaluate the BODC for each isotopomer. We used the same values of 
/?, r e , and 0 e as in our previous work, [3] and included all terms with i + j + k < 6. It should 
be noted that was erronously described[3] as having the value 2 a 0 2 . It actually has the 
value 2 A~ 2 . The rms errors for the fits to F 0 were 1.8 x 1CT 5 au and 7.7 x 1CT 5 au, resulting 
in rms errors in AV of less than 0.003 cm“l. The linear coefficients cf jk were determined by 
equally weighted least squares fits. Due to numerical errors, the fits did not quite have the 
desired symmetry, so the coefficients were symmetrized after the fitting. 


IV. GENERALIZING THE HANDY METHOD 


In this section we generalize the Handy method for the BODC to enable one to compute 
the non-adiabatic correction. We start first with the kinetic energy operator for the nuclei 
in an arbitrary space fixed frame: 



Q 



d 2 


(4) 
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where m a is the mass of nuclei a, and i runs over cartesian components. We can re-write 
this in the form 

( 5 ) 

2 ^{dX, 


m a 0X ta 1 

where f means complex conjugate acting to the left. Now we can compute the non-adiabatic 
correction functions as reported previously, [5] using the space fixed cartesians of Exp 4. In 
this case one obtains the correction 
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( 12 ) 

( 13 ) 


where W n is the Born-Oppenheimer energy of electronic state n, and ip n is the Born- 
Oppenheimer electronic wavefunction for electronic state n. When computing the deriva- 
tives, the origin used for the electronic wavefunctions is fixed. It should be noted that in 
general, the cO Q ,, , and C (0) are dependent on the geometry of the nuclei. 

It will be convenient to express the coefficients in Eq.G :.n terms of mass independent 
quantities. Now m a m a CO Q > is independent of mass, so that part is straightfoi waid. Foi 
the other quantities, we need a partial decomposition, so that 
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and 


1 1 


(15) 


C (o) = E 




i rn Q m a > 

The transformation from atomic cartesians to internal coordinates with total angular 
momentum eigenfunctions takes place in several steps. We now considei each step in turn. 


A. Mass factor transformation 

We next introduce the coordinate transformation 

X = xM, 


(16) 


where x is the matrix of cartesian internal coordinates, and M is a mat rix of mass factors. [15 
17] We assume that the final vector is for the center of mass, and M is invertible. Let M 
denote the inverse of M. Then since 
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where (3 indexes the internal coordinates, the kinetic energy operator becomes 
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B. Decouple Center of Mass 


We next perform a transformation which takes the electronic coordinates and the coordi- 
nates for the nuclear center of mass and forms new coordinates consisting of the electronic 
coordinates with their origin the nuclear center of mass, anc. the coordinates for the total 
center of mass. [5] This transformation zeros the coupling between the center of mass and 
all other coordinates, and requires the inclusion of the mass polarization t.erm[5] when com- 
puting matrix elements. Nevertheless, the numerical values of the matrix elements and the 
remaining correction functions are unchanged. 

C. Euler internal 

We next transform from the space fixed nuclear cartesians to body fixed coordinates 
defined by 

X = A t x so , (24) 

with A a rotation matrix[17] parameterized by the Euler angles ad 7, the superscript T 
denotes transpose, and x s0 the z embedding of Sutcliffe and Tennyson, [18] parameterized 
by the lengths n and r 2 , the angle x, and the embedding parameter a. The matrix A can 
be easily determined by the relation 

A = V(x so )V(x) 7 , (25) 

where V is a 3 x 3 orthogonal matrix constructed from the first two vectors of its argument. 
We take the first column of V to be the unit vector parallel to the first vector of its argument, 
the third column to be the unit vector obtained by Schmidt orthogonalizing the first vector 
from the second, and the second column the cross product of the other two columns. 

Then we introduce[19] 

= (26) 

j l 

with Oi and g m given in Table I. The coordinates / consist of the radial coordinates n and 
r 2 , the bending angle x, and the nuclear angular momentum operators R^. In the table, R ^ 
is the £ component in the nuclear body-fixed frame of reference of the total nuclear angular 
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momentum, with its sign reversed to satisfy normal commutation relations. Substituting 
Eq. 26 into Eq. 18 easily yields the desired kinetic energy operator: 

T=^yO,'G„.O r , (27) 

1 ll' 


with 
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(28) 


We must also remember to include the volume element obtained from the Jacobian deter- 
minant when computing matrix elements over vibrational-rotational functions. This volume 


element is not to be acted on by any of the O l when T is written as in Eq. 27. We can also 
write Eq. 27 in the more conventional form 


r = %Y.GivoA, 

~ U ' 


(29) 


As before, we obtain the noil-adiabatic correction 
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(33) 


D. Rotational invariance 

Finally, we fix the electronic positions with respect to the nuclear body-fixed axis. This 
results in the elimination of the nuclear angular momentum operator via the relation 

/?£ = , (33) 

where £ specifies the cartesian component in the nuclear body-fixed frame of reference, R 
is the nuclear angular momentum with sign reversed to satis :v normal commentation rela- 
tions, J is the total angular momentum with sign reversed to satisfy normal cominuntation 
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relations, and L is the electronic angular momentum. Since the electronic wavefunetion 
rotates with the molecule. Jc will have no effect on the electronic wavefunctions, thus ma- 
trix elements which were computed as derivatives with respect to the Euler angles ait now 
computed as matrix elements of L$, but the numerical values are unchanged. 


E. Wavefunetion Factorization 


So far we have presented results where the radial coordinates have volume elements r‘ 2 rfr,. 
It is convenient to introduce the factorization 4 ; ;; = ^/(ri r-i) and solve for T,, rathei than 
for T,,. This changes some of the coefficients in Eq. 30. Specifically ] is replaced with 


C m - 2 cS'Vn - 2C‘ l) /r. 2 + cll’/rf + dStA + 2 C^/nr, 


(35) 


di) . 


and then C t is replaced by 
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—12 / ' I' 


(36) 


F. Summary 

In summary, one computes m a m n >C^j a ', C^ a ,, using the appropriate modifications 
of Eqs. 7,8,9 on a grid. Then for a particular isotope, one forms c£\ and C (0) 

(Eq. 14,15), then applies the mass transformation of Eqs. 21 and 22, then the rotation and 
coordinate transformation of Eqs. 31 and 32, and finally the factorization transformation 
of Eqs. 35 and 36. One could attempt to fit m Q nvC^ a , etc directly, but since these 
quantities depend on the choice of origin, it is not clear what functional form to use. For 
this part of the calculation, we fit the C?, and C (0) calculated at 366 geometries with 
165 parameters to a functional form like the one described previously. [5] 


V. RESULTS 


Now that we have the icMRCI BODC, we can test how well the BODC computed using 
the CASSCF wavefunetion [5] is converged. In Table II we give our results obtained using 
the SCF BODC, the CASSCF BODC, and our new icMRCI BODC. All calculations used 
the cc-pVTZ basis. [20] These calculations use the un-optimized PES V 5Z as in Ref.5. To 
be consistent with our previous work. [5] we label the vibrational levels with local mode 


11 


quantum numbers. We see that we obtain quite similar results from all methods for bending- 
overtones. but much larger differences for stretching overtones. The corrections obtained 
from the CASSCF method are overestimates for the stretches, for the icMRCI results always 
move toward the SCF results. Nonetheless, the SCF results are not reliable for stretches. 

We now turn to our results using the new PES. In Table III, we give the errors for low 
lying levels of the various isotopomers from three different calculations: the results labeled 
ymass are taken f rom R e f. [3] and are obtained using the empirical mass dependent correction, 
the results labeled BODC use the PES V emp from Ref. [3] and the icMRCI BODC from the 
present work, and finally the results labeled Noad use the newly optimized PES V nemp and 
include the icMRCI BODC and the non-adiabatic corrections. Since the PES did not include 
the BODC in its optimization, we subtract the H 2 16 0 BODC from the BODC computed 
for the isotopomer of interest. In this table, the levels are labeled by normal mode quantum 
numbers. Comparing the results from BODC and Noad, we see that in all cases but one, the 
Noad results are an improvement over the BODC results, and in some cases significantly so. 
The exception is the IIO stretch in HTO, but both calculations do quite well for this level. 

We can see how close we are to our goal of accurate ab initio predictions of isotopic- 
shifts by comparing the errors to the average errors obtained for H 2 16 0, which was used to 
optimize the PES. The rms errors for the H 2 16 0 levels in Table III is 0.04 cm -1 . There are 
five errors larger than this in the table for the other isotopomers. Of the five, three are for 
the anti-symmetric stretch of D 2 16 0, T 2 16 0, and HO stretch of HT 16 0, and the largest is 
-0.10 cm -1 . These are probably not significantly larger than we expect. The remaining two 
are larger: 0.18 cm" 1 for the symmetric stretch of D 2 16 0 and 0.22 enr 1 for the symmetric 
stretch of T 2 16 0. Since these errors are larger and are for the same mode, we think that 
these are real, significant, discrepancies. We are not sure of the origin of these differences, 
but we suspect that it is due to a deficiency in our calculations. 

It is interesting to also compare our present best results to the previous ones using the 
empirical mass correction. In all but four cases, the empirical correction gives smaller cirois. 
In these four cases, the errors are all small, except for the HO stretch in HT 1 ’O, where 
ymass gj ve s —0.14 cm -1 and our best ab initio results are —0.10 cm -1 . In contrast, in the 
two cases mentioned above where the ab initio calculations gave large enors, the empiiical 
mass correction gives significantly smaller errors: 0.00 cm 1 and —0.11 cm h However, on 
the whole, the two results are very similar, which is a very encouraging result. 
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VI. CONCLUSIONS 


We have computed the isotopic shifts for the water molecule including both the diagonal 
and non-adiabatic corrections to the Born-Oppenheimer approximation. We find that both 
corrections are significant, and on the whole, the isotopic shifts are predicted quite well for 
low lying levels. The exceptions are the symmetric stretching fundamentals of D 2 16 0 and 
T 2 16 0, where the errors are about 0.2 cm -1 . 

There are four possible sources of the disagreement: (i) the Born-Oppenheimer PES is 
not accurate enough, i.c. cancellation of errors allow most of the isotopomeis to be fitted 
accurately, but. not all, (ii) the BODC is not as accurate as we think possible sources of 
error are the use of the cc-pVTZ basis and icMRCI wavefurctions, (iii) the non-adiabatic 
correction is not as accurate as we would like, cither due to the inability of the scaling to 
make up deficencies in the SCF/CIS method, the basis set is too small so that incoirect 
geometry dependence is obtained, or the fit is not good enough, or (iv) the experimental 
data is incorrect. Work is ongoing to test these possibilities. We are attempting to compute 
the best possible ab initio PES [2 1 ] as w'ell as generalizing our non-adiabatic code to more 
accurate wavcfunctions. Nonetheless, the future is very bright for accurate predictions of 
isotopic shifts. 
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TABLE I: Coordinate transformation 

matrix. 


d d 

dr i dr? 

d_ 

Ox 

R x 

Ry 

R z 

— sin ax 0 

- cos ax/n 

0 

(1 — a) cos ax/ir i 

0 

0 0 

0 

-sin(l - a)x/(*>i sinx) 

0 

— cos(l — a)x/(*n sin 

cos ax 0 

-sin ax/n 

0 

— a) sinax/^‘i 

0 

0 sin(l— a)x 

cos(l - a)x/ r 2 

0 

acos(l — a)x/ ir 2 

0 

0 0 

0 

— sin ax/ (in sinx) 

0 

cosax/(in sinx) 

0 cos(l — a)x 

- sin(l - a)x/n 

0 

-a sin(l — «)x/ 2 * r 2 

0 
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TABLE II: BODC for II 2 0 from various electronic wavefunctions (incm '). 


Ws 

SCF 

CASSCF 

icMRCI 

0(0,0) + 

0.00 

0.00 

0.00 

1(0,0) + 

-0.50 

-0.4G 

-0.45 

2(0,0)+ 

-1.00 

-0.94 

-0.91 

0(1,0)+ 

-0.05 

0.39 

0.28 

0(1,0)- 

0.13 

0.62 

0.47 

3(0,0)+ 

-1.52 

-1.46 

-1.40 

1(1,0)+ 

-0.49 

-0.01 

-0.11 

1(1,0)- 

-0.31 

0.24 

0.09 

4(0,0)+ 

-2.08 

-2.03 

-1.95 

2(1,0)+ 

-0.94 

-0.42 

-0.50 

2(0,1)- 

-0.74 

-0.16 

-0.29 

0(2,0)+ 

-0.05 

0.96 

0.67 

0(2,0)- 

0.09 

1.19 

0.83 

0(1,1)+ 

0.21 

1.22 

0.90 
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TABLE III 

: Errors and 

experimental energies 

for low-lying 

vibrational love 

ds of various isotopic 

substitutions (in cm 1 ). 





Source 

h 2 16 o 

hd 16 o 

d 2 16 o 

t 2 16 o 

ht ,6 o 




010 



ymass 

—0.03 

-0.04 

-0.03 

0.00 

-0.01 

BODC 


0.02 

0.00 

-0.04 

-0.04 

Noad 

0.05 

0.01 

0.00 

-0.02 

0.01 

exp. 

1594.78“ 

1403.48“ 

1178.38 6 

020 

995.33“ 

1332. 48 d 

yrnass 

0.00 

-0.01 

0.00 



BODC 


0.20 

0.06 



Noad 

0.05 

0.07 

0.00 



exp. 

3151.63“ 

2782.01“ 

2336.84“ 

100 



y mass 

0.01 

0.02 

0.00 

-0.11 

-0.04 

BODC 


0.51 

0.62 

0.39 

0.31 

Noad 

—0.04 

0.00 

0.18 

0.22 

0.03 

exp. 

3657.04“ 

2723.68“ 

2671.65“ 

001 

2237.15 7 

2299.779 

y mass 

-0.03 

-0.01 

-0.01 

-0.01 

-0.14 

BODC 


0.76 

0.58 

0.29 

-0.07 

Noad 

0.03 

0.02 

-0.07 

-0.08 

-0.10 

exp. 

3755.96“ 

3707.47“ 

2787.72“ 

2366. 60^ 

3716. 58-^ 


“Ref. 7 
6 Ref. 22 
“Ref. 23 
d Ref. 24 
e Ref. 25 
^Ref. 26 
9 Ref. 27 
A Ref. 28 
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